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1.  INTRODUCTION 


This  note  presents  the  solution  formulation  and  finite  element 
discretization  of  a  stress  wave  problem  with  discontinuous  data  in  two 
variational  schemes.  The  first  is  in  a  sense  a  generalized  Galerkin's 
approach  in  that  it  works  for  non-self  adjoint  problem  and  that  all  the  end 
conditions  are  made  to  be  natural  ones  and  hence  none  of  them  are  required  to 
be  satisfied  by  the  trial  functions.  This  unconstrained  variational  finite 
element  formulation  has  been  applied  to  initial/ boundary  value  problems  other 
than  wave  equations  previously . 1 » 2  More  recently,  an  adjoint  bilinear 
variational  principle  has  been  developed  for  initial  and  initial/boundary 
value  problems  which  requires  that  the  initial  conditions  oe  satisfied  exactly 
and  the  variations  of  the  adjoint  variable  be  set  to  vanish. ^  It  is 
consequently  a  constrained  variational  formulation.  This  note  compares  one 
formulation  and  numerical  results  with  those  of  the  other. 

First,  in  Section  2,  the  physical  problem  of  a  longitudinal  stress  wave 
in  an  elastic  rod  is  stated.  The  rod  is  fixed  at  one  end  and  free  at  the 
other  end.  The  discontinuity  data  arises  from  the  initial  linear  displace¬ 
ment,  corresponding  to  a  constant  stress,  due  to  a  force  applied  at  the  "free" 


lj.  J.  Wu,  "The  Initial  Boundary  Value  of  Gun  Dynamics  Solved  by  Finite 
Element  Unconstrained  Variational  Formulations,"  Innovative  Numerical 
Analysis  For  the  Applied  Engineering  Science,  R.  P.  Shaw,  et  al,  Editors, 
University  Press  of  Virginia,  Charlottesville,  pp.  733-741,  1980. 

2j.  J.  Wu,  "Solutions  to  Initial  Value  Problems  by  Use  of  Finite-Eleraents- 
Unconstralned  Variational  Formulations,"  1977  Journal  of  Sound  and  Vibration, 
53,  p.  341-356. 

■^C.  N.  Shen  and  J.  J.  Wu,  "A  New  Variational  Method  for  Initial  Value 
Problems,  Using  Piecewise  Herraite  Polynomial  Spline  Functions,"  presented  at 
the  1981  Army  Numerical  Analysis  6  Computers  Conference,  Huntsville,  AL, 
February  1981. 
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end.  This  force  suddenly  disappears  at  time  zero  causing  a  stress  discontinu¬ 
ity  at  the  free  end.  The  two  variational  formulations  for  the  stated  problem 
are  Introduced  in  Section  3.  finite  element  discretization  and  shape  func¬ 
tions  are  introduced  in  Section  4.  finally,  numerical  results  and  comparisons 
are  made  In  Section  3. 


2.  STATEMENT  OF  THE  PROBLEM 

The  problem  considered  here  is  that  of  a  longitudinal  stress  in  a  rod. 
The  differential  equation  can  be  written  as 

a2u  i  <i:u  0  s  *  s  x 

___  =  _  — -  ;  ( n 

3x‘  ot"  o  S  t  s  r 

with 

a*-  =  E/p  (2) 

where  u  ®  u(x,t)  is  the  axial  displacement 

x,t  are  the  coordinates  in  axial  direction  and  In  time,  respectively 
p  ,E  are  density  and  Young's  modulus,  respectively,  of  the  rod  materia) 
t  »  length  of  the  rod 
T  =»  some  finite  time  of  Interest 

For  the  boundary  conditions,  we  wilL  consider  a  rod  fixed  af  one  end  and 
not  restrained  at  the  other  end.  Hence 

u(0,t)  »  i 

3u  <  ” 

— ( i  ,t )  -  :> 

3x 

The  dynamics  of  the  problem  is  due  to  the  initial  conditions.  It  is 
assumed  that  the  rod  Is  stretched  to  a  linear  displacement  by  a  force  P  which 
vanishes  at  time  t  >  0  (see  Figure  1).  The  initial  velocity  of  the  rod  is 
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assumed  to  be  zero.  Thus 


P 

u(x,0)  -  —  x 
AE 

3u  O) 

— (x.O)  =  0 

at 


Figure  1.  Problem  Configuration  and  Applied  Load  at  Zero  Time 
(i.e. ,  P  -  0  for  t  >  0). 

It  Is  convenient  to  use  dimensionless  parameters.  Let 

u  =*  u/1  ,  x  -  xM  ,  t  =*  t/T  (5) 

Then,  Eq.  (1)  In  dimens  ionless  form  Is 


where 


9211  v,2 

3zu 

0  <  x  <  1 

9xz  = 

3tz  ' 

0  <  t  <;  1 

(b) 


(7) 


The  boundary  conditions  become 

u(0,t)  *■  0 


9u 

—  O.t) 
3x 


0 


(«) 
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and 


-  -  —  3u 

u(x,0)  -  Px  ,  —  (x  ,0)  -  0  (9) 

3t 

where 

P 

P  -  —  (10) 

AE 

is  the  force  in  dimensionless  form. 

The  stated  problem  in  dimensionless  form  are  Eqs.  (6),  (8),  and  (9)  with 
the  new  dimensionless  parameters  related  to  physical  counterparts  by  Eqs.  (5), 
(6),  and  (10).  To  simplify  writing,  we  shall  drop  the  bars  in  Eqs.  (6),  (8), 
and  (9),  and  rewrite  them  as 

0  <  x  <  1 

u"  -  b2u  *  0  ;  ( 6 ' ) 

0  <  t  <  1 

u(0,t)  -  0  ;  u'(l,t)  *  0  (S’) 

u(x,t)  -  Px  ;  u(x,0)  *  0  •  (9') 

where  a  prime  (')  indicates  differentiation  with  respect  to  x  and  a  dot  (•), 
with  respect  to  t. 

3.  TWO  VARIATIONAL  FORMULATIONS  OF  SOLUTIONS 
Consider  a  variational  problem 

6I0  *  0  (11a) 

with 

iQ  *  Iq(u,v)  -  /  J  (-u’v'+b2uv)dxdt  (lib) 

0  0 

where  u(x,t)  and  v(x,t)  are  said  to  be  adjoint  to  each  other.  It  is  a  simple 
matter  to  see  that  this  problem  is  an  indeterminate  one.  However,  the 
functional  of  Eq.  (lib)  can  be  modified  to  a  variational  problem  which  is 
equivalent  to  the  boundary/initial  problem  of  Eqs.  (6'),  (8'),  and  (9'). 
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Thus,  consider 


61  -  0  (12a) 

with 

,1  ,l  ,*• 

I  =*  I(u,v)  “  i  j  (~u’v  '+h“uv  )dxdt 
0  0 

,1 

+  ki!  u(0,t )v(0,t )dt 

0 

+  kob2/  [u(x,0)  -  u0(x) ]v(x ,  l)dx  +  b2  /  ui(x)v(x  ,0)dx  (12b) 

0  0 


We  shall,  take  the  first  variation  of  the  function  I(u,v)  of  Eq.  (12b)  in 
such  a  manner  that  6v  Is  completely  arbitrary  while  6u  is  set  to  zero 
Identically.  Hence,  by  means  of  integrations-by-parts ,  one  has 

,1  ,1  •' 

(6I)6u»0  3  JQ  JQ  (u"-b2u)6vdxdt 

.1 

-  u' ( 1 ,t)6v( 1 ,t)dt 

,1 

+  [u(0,t)  +  kiu(O.t) ]6v(0,t)dt 

+  b2J  (u(x , 1 )  +  k2[u(x,0)  -  uo(x) ] } 6v(x  ,  1  )dx 
0 

-  b  J q  [u(x,0)  -  ui(x) ) 6v(x,0)dx  =  0  (13) 

The  fact  that  5v(x,t)  is  completely  arbitrary  enables  us  to  conclude  from  Eq. 
(13)  that 

0  <  x  <  1 

u"  -  b2u  =  0  ; 

0  «  t  M 
u’(l.t)  =  0 

u'(0,t)  +  klu(0,t)  -  0  (14) 

u(x,l)  +  k2lu(x,0)  -  u0(x)  ]  =*  0 
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and 

u(x,o)  -  m (x )  -  o 

It  is  then  observed  that  the  initial/boundary  value  problem  defined  by  Eq. 

(14)  reduces  to  that  of  Eqs.  (6'),  ( 8 * ) ,  and  (9')  if  one  lets  kj  and  k.2  go  to 
infinity  (and  with  uG(x)  -  Px  and  u[(x)  -  0).  This  fact  suggests  that  the 
variational  problem  of  Eqs.  (12)  can  be  used  as  a  basis  of  a  finite  element 
discretization  for  the  approximate  solutions  to  the  original  initial/boundary 
problem.  It  should  be  noted  that  all  the  axiaiiary  conditions  in  Eqs.  (14) 
are  the  so  called  natural  boundary  conditions.  They  are  the  consequence  of 
the  variational  problem  -  just  like  the  differential  equation  itself.  For 
this  reason,  the  above  solution  is  referred  to  as  an  unconstrained  variational 
formulation. 

Another  approach  begins  from  Eq.  (lib).  With  6u  =  0  once  again,  one  has 

,1  -1  • 

6I0  -  /  (-U1 6v'+b‘-u5v)dxdt 

rl  -1 

+  /  u'(  1  ,t)<5v(  1  ,t)dt  -  j  u' (0  ,t )  6v(  0  ,t )dt 

0  0 

-  b2/  u(x ,  l)6v(x ,  1  )dx  +  b2  j  u(x ,0) 5v(x  ,0)dx  =»  0  (15) 

0  0 


with  the  constrained  conditions 


u(0,t)  -  0  ;  ul(l,t)  »  0  for  0  <  t  <  1 

(16) 

u(x,0)  “  u0(x)  ;  u(x,0)  =«  0  for  0  <  x  <  1 

It  was  shown  in  another  paper^  that  the  variations  of  the  adjoint  variable 


^C.  N.  Shen,  "Method  of  Solution  for  Variational  Principle  Using  Bicubic 
Hermite  Polynomial,"  presented  at  the  17th  Conference  of  Army  Mathematicians, 
West  Point,  NY,  June  1981. 
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must  he  constrained  as  follows 


£v(l,t)  *  0  ;  ov'(0,t)  =  0  for  -i  x  l  x  j 

(17; 

5v(x,l)  "  1  ;  cvix.l)  •=  0  for  0  >■  x  ^  i 


4.  FINITE  ELEMENT  DISCRETIZATION  AND  SslAFK  Fl  NCT  I ' ' N  S 

Through  non.limenslon.il  1  zat  ton  ,  the  region  o:  l:.t  -rest  it  wa  renal  ns  to 
be  a  unit  square:  0  x  x  <  1  and  Os  t  •*  1.  The  finite  element 
discretization  is  a  subdivision  of  this  unit  square  into  smaller  rectangles, 
the  elements.  A  typical  element  scheme  Is  shown  in  rirutv  2  where  a  typical 
(i,j)t*'1  element  Is  also  shown .  In  terms  of  me  element  car iahien  5q.  (  12  t)  is 


now  written  as 


51 


*•  ;i(t,j) 

f.i 


(IS) 


Variables  u(x,t)  and  v(x,t)  become  U(  j.  «  )  (  y  ,  -  }  and  i  ,  i  )  (  , '  )  respectively 
where  ,n  are  local  independent  variables  in  spatial  and  temporal  axis  also 
shown  in  Figure  2. 

Relations  between  global  and  local  coordi nates  are  >.:  ven  ;<  f  il  lows 

C  =  =  Kx  -  i  h  ! 


•  19) 


h  =  r,(J  )  =  Lt  -  j  t  1 

where  K  and  L  are  the  number  of  segments  In  x  and  t  directions,  respectively 
(see  Figure  2). 

Shape  functions  are  Introduced  as  follows.  Let 


u(i,j)(5,n)  =  aT(r,,o)U({  j) 


(20) 


where  a(£,n)  Is  the  shape  function  vector  and  U(t,j)  Is  the  discretized 
unknown  vector.  In  this  paper,  a(£,h)  is  selected  as  the  following. 
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With  the  conventions  as  stated  above,  the  meaning  of  the  unknowns  U)<(i,j)  in 

the  vector  U(i,j)  is  as  follows 

3u  3u  32u 

Ui  -  ,u(0 ,0)  ;  U2  -  —(0,0)  ;  U3  -  —(0,0)  ;  U4— — (0,0) 

35  3n 


3u 

3u 

-  —(0,0) 

;  u3  -  — (0, 

H 

3n 

3  u 

3u 

— (1.0)  ; 

U7  -  —  d.O) 

31 

3r. 

3u  3u  3  u 

u9  -  u(o ,  1)  ;  U10  -  —(0,1)  ;  u  1 1  -  “(0,1)  ;  U12  -  ~ “(o.D 

35  3n 

U13-u(l,l)  ;  U14  -  ^(  1,1)  ;  Ui5  -  ~( 1,1)  :  Ul6  -  — (1.1) 

for  each  element  (i,j). 

For  the  unconstrained  formulation,  Eq.  (20)  is  used  in  Eq.  (12).  The 
result  is 


l  I  «y(1,j)T{-  7  A  +  b2  £  B}y(ltj)  +  [  5y(1,j)T(--)  c  u(ltj) 
1-1  j-l  l  k  j-i  L 


L  '  K 

L  k2b2  K  b^k2 

+  l  «y(i,i)T(--“)p  U(ifl)  -  l  5V(i,L)T(— ■ —  )F(1) 
i-1  K  1-1  K 


where 


K  b 2 

l  6-v(i,i)  (;~>9(i) 

i-i  K 


6  ■  / Q  JQ  ?,5§T,Cd^dri  ;  §  ”  JQ  /0  f,n?T,nd^dr| 


C  -  f1  a(0,n)aT(0,n)dn  ;  D  -  /  a(l,l)aT(l,0)dn 

0  ‘  '  ~0" 


F(t)  -  /  , l>dt  ;  ?(t)  -  /  U^i)a(^,0)d5 
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The  expression  of  F(i)  and  G(i)  .can  be  further  reduced  Into  a  form  more 


readily  computed.  Write 

u^^U)  -  aTU.O)^1) 


16 

I  ak(C.0)Uoic(1) 
k-1 


16 

U1U)(?)  -  aT  n(C.0)Uo<1>  -  l  ak,n<t.°)uok(1) 

k-I 


(26) 


Since 

akU.O)  -  bi(Obj(0) 
ak,nU,0)  “  biCOb’jCO) 

and 

b i(0)  -  1  ,  bj (0)  =  0  for  j  -  2,3,4 
b’2(0)  -  1  ,  b'j(O)  “0  for  j  -  1,3,4 
From  Table  I,  one  then  observes  that 

ak(5,0)  “  0  for  all  k  except  k  •  1,2, 5, 6 
ak,n^»°)  “  0  for  all  k  except  k  -  3, 4, 7, 8 
Hence,  In  Eq.  (26),  only  U0i(1>,  U02(1).  U()5(1)»  and  u06(1)  are  u8ed  ln 
expressing  Uq^^  and  only  Uq3^^  ,  G()4^^>  and  ^08^^  are  used  in 

expressing  ui'^(5)»  Thus  we  shall  write 


with 


F(i)  "  J*  a«,l)4TU,0)d&  UqU)  -  F  U,^1) 

*v  0  *  "  (27) 

9(1)  “  /q  a(C.0)aT>n(C,0)dC  U^1)  -  G  Uq^) 

F  -  J1  a(C,0)aT(c,0)dC 

0  ' 

(28) 

9  “  /q  aU,0)aT>nU,0)dS 


11 


The  way  to  aet  up  is  that  first  set  all  UQk^  to  zero  for  all 

k  -  1,2,. ..,16.  That  set  U0k^  for  k  -  1 , 2 , 3,4 , 5  ,6 , 7  and  8  as  follows. 
Uoi^1)  •  Uq^Ho)  ;  Uq2^1^  "  u0,£^(°)  »  Uo3^  “  uj^Ho)  ; 

U04(1)  “  ui,c(0)  ;  U05  “  vi0(l)(l)  ;  U06(l)  “  uO,t(l)<l)  ; 

Uo7^}  “  u1(i>(  1)  ;  Uo8^^  “  ul,4^^(l)  (29) 

With  vectors  U0(f),  F,  and  G  completely  defined  above,  Eq  .  (24)  can  he 

rewritten  as 


K  L  K  L  L  kt 

l  l  6V(i,J)T{-  -  A  +  h2  -  §}  U(lfJ)  +  l  6V(  1  ,j  )T(” )C  U(  1  ,j  ) 
1-1  j-1  L  K  j-1  L 


b2k2 


b2k2 


(30) 
(i) 


+  l  6-v(lfJ)T(---)p  u(lfl)  -  I  [6y(i>L)T(---)F  -  sy(i>1)T(--)r,]y0 

i- 1  K  i « 1  K  it 

Now  Eq.  (30)  is  readily  assembled  into  a  global  matrix  equation  in  a  standard 
manner. 

6VT  K  U  -  SVT  P  (31) 


or 

KU  -  P  (31) 

due  to  the  fact  6V  is  completely  arbitrary.  ThuB  Eq.  (31)  is  solved  for  U. 


5.  NUMERICAL  RESULTS  AND  COMPARISONS 


Some  preliminary  results  of  computation  are  presented  here.  de  shall  set 


b  -  1  in  the  differential  equation  (6')  for  simplicity.  Thus, 

,  P  %  2  *2 

b2  -  -(-)  -  -T-r  -  1  (32) 

E  T  azT2 


T  -  -  (33) 


The  exact  solution  for  t  -  0,  0.2T,  and  0.4T  are  given  in  Figures  3  and  4. 
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First,  the  results  from  thg  unconstrained  variational  formulation.  Using 
a  grid  scheme  of  Sxl,  the  numerical  results  for  displacements  and  axial 
stresses  are  tabulated  in  Tables  II  and  III  where  the  exact  solutions  are  also 
given  for  comparison.  The  graphic  comparisons  are  shown  in  Figures  5  and  6 
where  the  calculated  solutions  are  indicated  by  crosses  (x)  and  the  exact 
solution  is  plotted  in  solid  lines.  It  is  clear  in  these  figures  that  the 
computed  results  generally  agree  with  the  exact  analytical  solution.  As  a 
further  evidence  of  convergence,  a  finer  grid  scheme  of  10x1  is  taken  and  the 
improved  solution  is  shown  in  Figures  7  and  8. 

TABLE  II.  SOLUTIONS  TO  THE  STRESS  WAVE  PROBLEM  USING 
UNCONSTRAINED  VARIATIONAL  FORMULATION 
l 

t  -  0.2T  -  0.2(-),  b  -  1.0;  Grid:  5x1 
a 


u 

3u/  8x 

X 

Computed 

Exact 

Computed 

Exact 

0 

0.000 

0.0 

0.994 

1.0 

0.2 

0.199 

0.2 

0.989 

1.0 

0.4 

0.399 

0.4 

0.965 

1.0 

0.6 

0.598 

0.6 

0.896 

1.0 

0.8 

0.789 

0.8* 

0.550 

0.0* 

1.0 

0.806 

0.8 

0.403 

0.0 

Point  of  discontinuity 


TABLE  III.  SOLUTIONS  TO  THE  STRESS  WAVE  PROBLEM  USING  1 

UNCONSTRAINED  VARIATIONAL  FORMULATION  | 


t  -  0.4T  -  0 . 4 ( - )  ,  n  °  1.-  ;  GrU:  5x1 
a 


u 

3u, 

'3x 

X 

Computed 

Exact 

Computed 

Exact 

0 

0.000 

0.0 

0.988 

1.0 

0.2 

0.200 

0‘.  2 

0.976 

1.0 

0.4 

0.399 

0.4 

0.927 

1.0 

« 

t 

0.6 

0.594 

0.6 

0.714 

0.0* 

t 

0.8 

0.698 

0.6 

0  •  4  b  7 

0.0 

1.0 

0.791 

0.6 

0.151 

0.0 

*Point  ol 

discontinuity. 
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Figure  7.  Displacement  Solutions  by  the  Unconstrained  Variational 
Formulations  (t  =*  0.1T)  and  Comparison  with  Exact 
Solutions.  Grid:  10x1. 


] 
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Numerical  results  from  the  constrained  formulation  follow  the  general 
trend  as  the  unconstrained  one  as  Indicated  in  Figures  9  and  10,  as  well  as 
the  tabulated  comparison  with  the  exact  solution  in  Tables  IV  and  V. 

TABLE  IV.  SOLUTIONS  TO  THE  STRESS  WAVE  PROBLEM  USING 
CONSTRAINED  VARIATIONAL  FORMULATION 

t  -  0.1T,  b  -  1 ;  Grid:  5x1 


U 

9u/  3x 

X 

Computed 

Exact 

Computed 

Exact 

0 

0.0 

0.0 

0.986 

1.0 

0.2  . 

0.200 

o 

• 

ro 

0.984 

1.0 

0.4 

0.399 

0.4 

0.944 

1.0 

0.6 

0.596 

0.6 

0.784 

1.0 

0.8 

0.797 

0.8* 

-  0.049 

0.0* 

1.0 

0.874 

0.8 

0.0 

0.0 

*Point  of  discontinuity. 


I  I 


21 


I 


TABLE  V.  SOLUTIONS  TO  THE  STRESS  WAVE  PROBLEM  USING 
CONSTRAINED  VARIATIONAL  FORMULATION 

t  -  0.1T,  b  -  1.0;  Grid:  5x1 


'  u 

3u/  3x 

X 

Computed 

Exact 

Computed 

Exact 

0 

0. 

0.0 

1.000 

1.0 

0.1 

0.100 

0.1 

1.000 

1.0 

0.2 

0.200 

0.2 

1.000 

1.0 

0.2 

0.300 

0.3 

1.000 

1.0 

0.4 

0.400 

0.4 

0.999 

1.6 

0.5 

0.500 

0.5 

0.997 

1.0 

0.6 

0.600 

0.6 

0.986 

1.0 

0.7 

0.700 

0.7 

0.945 

1.0 

0.8 

0.798 

0.8 

0.787 

1.0 

0.9 

0.898 

0.9* 

-  0.038 

0.0* 

1.0 

0.936 

0.9 

0.0 

0.0 

* Point  or  discontinuity. 
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Figure  10.  Stress  Solution  by  the  Constrained  Variational 

Formulations  (t  =  0.2T)  and  Comparison  with  Fxact 
Solutions.  Grid:  5x1. 
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